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Abstract. We study a stochastic model of electricity production and consumption where appliances are 
adaptive and adjust their consumption to the available production, by delaying their demand and possibly us- 
ing batteries. The model incorporates production volatility due to rencwables, ramp-up and ramp-down time, 
uncertainty about actual demand versus planned production, delayed and evaporated demand due to adaptation 
to insufficient supply. We study whether threshold policies stabilize the system. The proofs use Markov chain 
theory on general state space. 



Key words. Dynamical Systems, Smart Grids, Elastic Demand, Macroscopic Model, Stability 



AMS subject classifications. 60J05, 93E15 

1. Introduction. Recent results on modelling the future electricity markets [TT] suggest 
that they may lead to highly undesirable equilibria for consumers, producers or both. A central 
reason for such an outcome might be the combination of volatility in supply and demand, the 
delays required for any unplanned capacity increase, and the inflexibility of demand which leads 
to high dis-utility costs of blackouts. Further, the use of renewable energy sources such as wind 
and solar increases volatility and worsens these effects [5]. 

Flexible load is advocated in [5] as a mechanism to reduce ramp-up requirements and adapt 
to the volatility of electricity supply that is typical of renewable sources. A deployments report 
[3] shows the feasibility of delaying air conditioners using signals from the distributor. Adaptive 
appliances combined with simple, distributed adaptation algorithms were advocated in [H [TJ; 
they are assumed to reduce, or delay, their demand when the grid is not able to satisfy them. 
Some examples might be: e-cars, which may have some flexibility regarding the time and the 
rate at which their batteries can be loaded; heating systems or air conditioners, which can delay 
their demand if instructed to; hybrid appliances which use alternative sources in replacement for 
the energy that the grid cannot supply. If the alternative energy source is a battery, then it will 
need to be replenished at a later point in time, which will eventually lead to later demand. 

The presence of adaptive appliances may help address the volatility of renewable energy 
supply, however, backloggcd demand is likely to be merely delayed, rather than canceled; this 
introduces a feedback loop into the global system of consumers and producers. Potentially, one 
might increase the backlogged demand to a point where future demand becomes excessive. In 
other words, one key question is whether it is possible to stabilize the system. This is the question 
we address in this paper. 

To address this fundamental question, we consider a macroscopic model, inspired by the 
model in [5]. We assume that electricity supply follows a two-step allocation process: first, in 
a forecast step (day ahead market) demand and renewable supply are forecast, and the total 
supply is planned; second, (real time market) the actual, volatile demand and renewable supply 
are matched as possible. We assume, as in pQ, that the rate at which supply can be increased 
in the real time step is subject to ramp-up and ramp-down constraints. Indeed, it is shown in 
this reference that it is an essential feature of the real time market. We modify the model in [5] 
and assume that the whole demand is adaptive. While this is clearly an exaggerate assumption, 
we do it for simplicity and as a first step, leaving the combination of adaptive and non-adaptive 
demand to a later research. We are interested in simple, distributed algorithms, as suggested 
in [4], therefore, we assume that the suppliers cannot directly observe the backlogged demand; 
in contrast, they see only the effective instantaneous demand; at any point in time where the 
supply cannot match the effective demand, the backloggcd demand increases. 
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Our model is macroscopic, so we do not model in detail the mechanism by which appliances 
adapt to the available capacity; several possible directions for achieving this are described in [3] . 
However we do consider two essential parameters of the adaptation process. First, the delay 1/A 
is the average delay after which frustrated demand is expressed again. Second, the evaporation [i 
is the fraction of backlogged demand that disappears and will not be resubmitted per time unit. 
The inverse delay A is clearly positive; in contrast, as discussed in Section [5] and in Appendix lAl 
it is reasonable to assume that some adaptive appliances naturally lead to a positive evaporation 
(this is the case for a simple model of heating systems), but it is not excluded that inefficiencies 
in some appliances lead to negative evaporation. 

Within these modelling assumptions, the electricity suppliers are confronted with a schedul- 
ing issue: how much capacity should be bought in the real time market to match the adaptive 
demand. The effect of adaptation is to increase the latent demand, due to backlogged demand 
returning into the system. This is the mechanism by which the system might become unstable. 
We consider a threshold based mechanism as in [8]. It consists in targeting some fixed supply 
reserve at any point in time; the target reserve might not be met, due to volatility of renewable 
supply and of demand, and due to the ramp-up and ramp-down constraints. 

Our contribution is to show that if evaporation is positive, then any such threshold policy 
does stabilize the system. In contrast, if evaporation is negative, then there exists no threshold 
policy which stabilizes the system. The case where evaporation is exactly equal to remains 
unsolved. 

Our results suggest that evaporation plays a central role. Simple adaptation mechanisms 
as described in this paper might work if evaporation is positive (as one may perhaps generally 
expect), but will not work if evaporation is negative, i.e. if the fact that demand is backlogged 
implies that a higher fraction of demand returns into the system. This suggests that future 
research be done in order to gain a deeper understanding of evaporation, whether it can truly 
be assumed to be positive, and if not, how to control it. 

We use discrete time, for tractability. We use the theory of Markov chains on general state 
spaces in [S] . In Section [3] we describe the assumptions and the model, and relate our model to 
prior work. In Section [3] we study the stability of the system under threshold policies. Proofs 
are given in Section 2] The appendix discusses whether evaporation can be positive or negative 
when the appliance is a heating system, including heat pumps. 

2. Model and Assumptions. 

2.1. Assumptions and Notation. We use a discrete model, where f e I represents the 
time elapsed since the beginning of this day. The time unit represents the time scale at which 
scheduling decisions are done, and is of the order of the second. 

The supply is made of two parts: the planned supply G$ (t), forecast in the day- ahead market, 
and the actual supply G a (t), which may differ, due to two causes. First, the forecasted supply 
may not be met, due to fluctuations, for example in wind and sunshine. Second, the suppliers 
attempt to match the demand by adding (or subtracting) some supply, bought in the real time 
market. We assume that this latter term is limited by the ramp-up and ramp-down constraints. 
We model the actual supply as 

G a (t) = G(t-l) + G f (t)+M(t) (2.1) 

where M(t) is the random deviation from the planned supply due to renewables, G(t — 1) is 
the supply decision in the real time market. We view G a (t) as deterministic and given(it was 
computed yesterday in the day-ahead market), M(t) as an exogenous stochastic process, imposed 
by nature, and G(t — 1) as a control variable. 
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We call D a (t) the "natural" demand. It is the total electricity demand that would exist if 
the supply would be sufficient. In addition, there is, at every time t, the backlogged demand B(t), 
which results from adaptation: B(t) is the demand that is expressed at time t due to a previous 
demand being backlogged. The total effective demand, or expressed demand, is 

E a (t) = D a (t) + B{t) (2.2) 

We model the effect of demand adaptation as follows. 

B(t) = \Z(t) (2.3) 
Z(t + 1 ) = Z(t) - B(t) - fiZ(t) + F(t) (2.4) 
F(t) = [E a (t)-G a (t)} + (2.5) 

In the above equations, F(t) is the frustrated demand, i.e. that is denied satisfaction at 
time t. Eq. (|2.5|) expresses that, through adaptation, the demand that is served is equal to 
the minimum of the actual demand and the supply. The variable Z(t) is the latent backlogged 
demand; it is the demand that was delayed, and might later be expressed. It is incremented by 
the frustrated demand. 

The frustrated demand is expressed with some delay; we model this with Eq. (|2.3[) . where 
A -1 is the average delay, in time slots. 

The evolution of latent backlogged demand is expressed by Eq. (|2.4[) . The expressed demand 
B(t) is removed from the backlog (some of it may return to the backlog, by means of Eq. (|2.5|) 
at some later time). The remaining backlog may evaporate at a rate /i, which captures the effect 
on total demand of delaying some demand. Delaying a demand may indeed result in a decreased 
backlogged demand, in which case the evaporation factor is positive. For example this occurs if 
we delay heating in a building with a heating system that has a constant energy efficiency (as we 
discuss in Appendix IA.1[) ; such a heating system will request more energy in the future, but the 
integral of the energy consumed over time is less whenever some heating requests are delayed. In 
this case, positive evaporation comes at the expense of a (hopefully slight) decrease in comfort 
(measured by room temperature). In other cases, though, we may not exclude that evaporation 
be negative. This may occur for example with heat pumps, as discussed in Appendix IA. 21 

As in [5], we assume that the natural demand can be forecast with some error, so that 

D a {t) = D{t) + D f {t) (2.6) 

where the forecasted demand D> (t) is deterministic and D(t), the deviation from the forecast, 
is modelled as an exogenous stochastic process. We assume that the day ahead forecast is done 
with some fixed safety margin ro, so that 

G f {t) = Df{t)+r a (2.7) 

2.2. The Stochastic Model. We model the fluctuations in demand D{t) and renewable 
supply M(t) as stochastic processes such that their difference M(t) — D(t) is an ARIMA(0, 1, 0) 
Gaussian process, i.e. 

{M(t + 1 ) - D(t + I)) - (M(i) - D{t)) = N(t + I) (2.8) 

where N(t) is white Gaussian noise, with variance a 2 . This is the discrete time equivalent of 
Brownian motion, as in [3]. 

Let R{t) be the reserve, i.e the difference between the actual production and the expressed 
demand, defined by 



R{t) = G a {t)~E a {t) 



(2.9) 
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and let H (t) be the increment in supply bought on the real time market, i.e. 

H(t) = G(t) - G(t - 1) (2.10) 

Putting all the above equations together, we obtain the system equations: 

R(t + 1) = R(t) + H{t) + N(t + 1) - A (Z(t + 1) - Z{t)) , (2.11) 
Z(t + 1) = (1 - A - fi)Z(t) + l {R{t)<0} \R(t)\, (2.12) 

Thus we can describe our system by a two-dimensional stochastic process X(t) = (R(t),Z(t)), 
with t £ IN. 

The sequence H(t) is the control sequence. It satisfies the ramp-up and ramp-down con- 
straints: 

-e < H(t) < c, 

where £ and £ are some positive constants. 

We assume a simple, threshold based control, which attempts to make the reserve equal to 
some threshold value r* > (; therefore 

H(t) = max (min (£ r* - R(t)) , -£) . (2.13) 

In summary, we have as model the stochastic sequence X = (X(t))t£N defined by 

R(t + 1) =R(t) - \l {m<0} \R(t)\ + A(A + n)Z(t) 

+ (C A (r* - i?(t))) V (-0 + iV(i + 1), (2.14) 
Z(t + 1) =l {fl(t)<0} |i?(<)| + (1 - A - n)Z(t\ (2.15) 

where is an iid white Gaussian noise sequence. Note that X is a Markov chain on the state 
space S = lx K+. 

2.3. Related Models. Let us now discuss some similar models which have been considered 
in the literature. 

In [5], the authors consider the so-called Linear State Space model (LSS), which introduces 
an n-dimcnsional stochastic process X = {Xk}k, with Xk 6 H™. For matrices F £ E, nx ™, 
G £ E, nxp , and for a sequence of i.i.d. random variables of finite mean taking values in W, the 
process evolves as 

X k = Flu + GW k , k > 1. (2.16) 

Our model (|2.14l) - (|2.15|) is in fact a superposition of three such LSS models, depending on the 
current state of the Markov chain. The challenge of showing that our model is stable comes from 
the fact that in the part of the state space in which R(t) < 0, the corresponding LSS does not 
satisfy the stability condition (LSS5) of [9] (which requires that the eigenvalues of F be contained 
in the open unit disk of C). 

In [TU] a slightly different model for capturing elasticity of demand is proposed. More specif- 
ically, the authors consider a scenario in which a deterministically bounded amount of demand 
arrives at each time step, while the supplier decides whether to buy an additional amount of 
energy from an external source at a certain cost. Unsatisfied demand is backlogged. A threshold 
policy is analyzed and found to be stable (the size of the backlog is found to be deterministically 
bounded) and optimal. Pricing decisions are also explored. The main differences with the present 
work are the following: 
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• Additional parameters which model delay and loss of backlogged demand (i.e. A and fi) 
enrich our model's expressivity. 

• We consider potentially unbounded demand, modeled as a 0-mean Gaussian random 
variable, which makes proving stability more challenging. 

• No results on pricing and cost-optimality are included in the present work. 

The continuous time model used in [5] seeks to capture the presence of two types of energy 
sources, primary and ancillary, the latter being less desirable (i.e. more costly) than the former, 
both subject to (different) ramp-up constraints. A threshold policy is again discussed in the 
context of rigid demand, which is simply dropped if not enough energy is available. The analyzed 
Markov chain is a two-dimensional process having on the first coordinate the quantity of energy 
used from the ancillary source in order to satisfy as much demand as possible, and on the second 
coordinate the reserve (i.e. energy surplus). 

3. System Stability under a Threshold Policy. Let us study how the system stability 
of (|2.14[) - (|2.15[) depends on the evaporation parameter fi. 
Define the following four domains: 

D l = (-oo,0) xR + , D 2 = [0,r*-C) x R+, D 3 = [r* -£r* +£) x R + , JJ 4 = [r*+£,oo) x R+. 

Then the process (|2.14p - (|2.15p rewrites in matrix form: 

4 

X(t + 1) = Hx(t)eD i} AiX(t) + t {X(t)£D 1 UD 2 }(o + l{X(t)eD 3 } r O _ l{x(i)eL> 4 }£o + N o(t + 1), 

i=l 

(3-1) 

where 



X(t) 



R(t) 
Z(t) 



,N (t) 



N(t) 




and 



A 1 = 



1 + A A(A + fi) 
-1 1 - A - fi 



1 A(A + ^) 
1 - A- fi 



,A, = 



A(A + fi) 
1 - A - fi 



The main reason for which the analysis of system stability is challenging is the fact that both 

A\ and Ai = A^ admit 1 as an eigenvalue. 

The first result of this section can be stated in the form of the following 

Theorem 3.1. If fi> 0, the Markov chain \2.1J$2/Hfy is positive Harris and ergodic. For 

any initial distribution p, the chain converges to its unique invariant probability measure in total 

variation norm, i.e. 



P (dx)I> n (x, .)-7T(.) 



The proof uses the theory of general state space Markov chains. The following Lemmas are 
instrumental in proving the result. The proofs are found in Section [4] 

Lemma 3.2. J/1 — A — fi < 1, then there exists a measure tp such that the Markov chain \2.1J\ 
\2.15\) is ip- irreducible. 

A set C is said to be i^p-small for some non-trivial measure vt and a positive integer T, if 
for all x £ C, the probability of reaching any measurable B is V T (x, B) > vt{B). 
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Furthermore, a set C is petite if there exists a distribution h on the positive integers and a 
non-trivial measure u^, such that for any x G C, and for any Borcl set £?, the transition kernel 
of the sampled chain has the following property: 



K h {x, B) := J2 h{t)F\x, B) > v h {B). 



A z^T-small set is implicitly 5-r-petite, where St is the Dirac distribution. 

Lemma 3.3. For any set C = Jx [0,6], with J a finite closed interval and b > 0, there exists 
Tq > and non-trivial measures (vt)t>t such that C is VT-small for all T > Tq. 

Two direct consequences follow: 

Corollary 3.4. The Markov chain \2. 14\t2.lS\) is aperiodic. 
Corollary 3.5. Any compact subset of the state space is petite. 

In our proof of Theorem 13.11 we exhibit a Lyapunov function for the system, which has 
negative drift outside of a compact set surrounding the origin. The existence of such a function 
along with the previous results, enable us to apply Theorem 13.0.1 of [9], by which we conclude. 

A second result is stated below: 

Theorem 3.6. If \i < 0, the Markov chain \2.l 4^2. 15\) is non-positive. 
In the proof, we find a function which has finite increments and which is such that outside 
a certain level set its drift is negative. By Theorem 11.5.2 from [5] we reach the conclusion. 
We sum up these results in the following 

Corollary 3.7. If a positive fraction of latent jobs disappear during each time slot (fi > 0J, 
then any threshold policy stabilizes the system. On the other hand, if delaying any job results 
in an increase of its requirement/workload by a positive fraction (fi < 0), then there exists no 
threshold policy that stabilizes the system. 

4. Proofs. Let 7 := 1 — A — fj, and note that 7 < 1 whenever [i > 0. 

4.1. Proof of Lemma 13.21 Fix some finite closed interval / and some a > 0, and consider 
the set E = / x [0, a]. 

Consider the measure ifiE defined as follows: for any Borel set A, 

ip E (A):=v(AnE), (4.1) 

where v denotes the Lebesgue measure on R 2 . Take a measurable set B C E. We denote the 
support sets of B on the two dimensions by 

Bi :={ret: (r,z) G B}, 
B 2 := {z G R+ : (r, z) G B}. 

We further denote the cross-section of B at some z G R+ by 

tt [z] B := {r G R : (r, z) G B}. 

Consider any initial point x = (r,z). We wish to show that whenever lpe{B) is strictly 
positive, there exists a finite T such that the probability of hitting B in T steps starting from x 
is also strictly positive. Then, by Proposition 4.2.1 (ii) from [5J, we can conclude. 

Let us thus consider Lp E [B) > 0. For any < e < 1, there exists a 5 > 0, such that the 
Lebesgue measure of the set B s := {(r, z) G B n E : z > 6} — B n E n (R x [5, 00)) is strictly 
positive and, furthermore, 

( = ^m = HB^_ 

<p E {B) ^e(B)' V ■ ; 
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For every < a < 1, we define the sets B (a) := {(r, z) : (r,az) G -B }. Clearly, 

K *V)) = ^ = ^, (4-3) 
a a 

and, furthermore, 

K[ z ]B s (a/3) = Tr [0z] B s (a), V < a, p < 1. (4.4) 
We consider an initial point x = (r, z) G 

In the following we consider a specific trajectory from x to B (more specifically to B s ) and 
we assess the probability of this trajectory. Let thus (R(0), Z(0)) = (r, z) and denote 

r(l) = (l + A)r + A(l- 7 )z + C, 
z(l) = — r + yz. 

We have that J2(l) = r(l) + AT(l) and Z(2) = 7 z(l) - r(l) - iV(l). 

We wish to select a large enough T such that the event u>\ = {N(l) < — r(l)} includes 
the event w 2 = {Z{2) G B2(7 T )}. We have that lu 2 C {72(1) - sup^( 7 T ) < r(l) + 2V(1) < 
7 z(l) — inf i?2(7 T )}- Since inf B^-y 7 ) > -4jr, it is sufficient to consider some T such that 

<* 



We namely pick a T such that T > Tq(B, e, x) A 0, where 



To(B,e,x) 



1 1 5 

■ loe 



log 7 z(l) 



(4.5) 



For such T, if we consider iV(l) such that u>2 occurs, then necessarily uj\ occurs, and R(l) < 0. 
In other words, at the beginning of time step 2 the process is still in region D\. 
For t = 2, . . . , T + 1, we consider N(t) such that R(t) G £(£), where 

L(e) = (r*-t,r*+£] 

for a fixed £ > £ A £. Since for such t the process evolves in domains D3 and D4, the evolution of 
Z(t) will be purely deterministic and can be written explicitly after 2 < t < T+ 2 steps as = 
7*" 2 Z(2). By definition, since Z(2) G B%(j T ), we have that Z(T + 2) G B|. We subsequently 
consider 7V(T + 2) such that (R(T+2), Z(T+2)) G £ 5 . Equivalently, (R(T + 2), Z(2)) G S 5 ( 7 T ). 

Let us now compute a lower bound of the probability of this trajectory. For some set Act 
and some real number b G R, we denote —A := {—a :aei} and ^4 + t>:={a + t>:aG A}. 

We can now characterize the probability of hitting B in T + 2 steps following the trajectory 
that we just described: 

V T+2 (x,B)> T(x,dxi) F(x u dx 2 ) ¥(x 2 ,dx 3 ) 

■/( 7 ^(l)-S 2 5 (7 T ))x{ Z (l)} JL^xB^T) Jl^xB^t-i) 



[ P(x T ,dxT+i)V{x T +uB s 



)xB 2 a (7) 

Denote the density of a normal variable of mean and standard deviation a by A/". 

The probability that N(l) is such that Z(2) G B^i^) (i.e., event 0J2) can be written as 

P{w 2 }= / JV(dm). 

7 7Z (l)-r(l)-B|( 7 ^) 
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We can lower bound the probability that N(2) is such that R(2) G L(£) as 

/.r*+«-C-(l+A)r-A(l-<y)*(l) 

p 2 (T,B,e,x):= inf / Af(dn 2 ) > 0, 

re 7 z(l)-B^(7 r ) 7 r »-2C-(l+A)r-A(l-7)z(l) 

since the infimum is taken over a bounded set. 

Similarly, after 2 < t < T+ 1 steps, the process will remain in D3 U (D4 n {(r, z) : r < r* +^}) 
with probability at least 

p t (T,B,e):= inf P e L{1) \ R(t - 1) = r, Z(t - 1) = z) . 

Note that R(t) has a different expression, depending on whether r < r* + £, or r > r* + ^. We 
introduce the following notation: 

f L(£)-r*-A(l- 7 )z, ifr <r*+£, 

1 ' ' \ L{1) - (r - £) - A(l - i)z, if r>r*+£. 

Then, we can rewrite 

Pt (T,B,e)= inf / N{dn t ) > 0. 

rGiW Jd(t,z) 
ze-B|( 7 T - f + 3 ) 

The strict positiveness of the lower bound follows from the fact that L(£) and B are bounded. 

In the last step we want to ensure that R(T + 2) belongs to the corresponding cross-section 
7T[ 7 T'(7 Z (i)-r(i)-n 1 )]S' 5 , which by dHU) is identical to B s {l T ). 

Let us denote 



(rznA-[ 7r [7-(i)-(i)-«i] B5 (7 T ) - r* - A(l - 7)*, i 



A. 

Then, we can finally write: 



if r < r* + £, 
if r > r* + £. 



P T+2 (x, 5) > ( TTpt] / AT(dn x ) inf / 

\t=2 / " / 7^(l)-''(l)--B2 i (7 T ) r£i W JD (r,z, ni 



N{dn T+2 ). 

(r.z.nx) 

7 T - 1 (7^(l)-r(l)-n 1 ) 

Since i(^) is bounded, there exists a large enough g such that 

inf / N(dn T+2 ) > / Af(dn T+2 ). 

reh(i) jD (r,x, ni ) Jniy. { l)-r(l)- ni ]Bt(l T )-q-Hl-y)z 

Z =7 T - 1 (72(l)-r(l)-n 1 ) 

We make the following coordinate change: 
v := 72(1) — r(l) — n\ 

w := n T+2 - q - A(l - 7)7 T_1 (7z(l) - r(l) - ux). 

Then, we can write: 



\t=2 J JJ B s h T ) 



w) dv dw, 
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where 

hx(v,w) = — exp |_-L[( 7Z (i) - r (i) - v f + (w - A(l - 7 ) 7 T - 1 v - qf 

Now since B s (j T ) is a bounded set, there exists r] = rj(T, B, e, x) > such that /^(v , w) > r/, for 
all (v,w) G B' 5 (7 T ). We can then conclude by (|4.3[) : 

P T+2 (x-, S) > ^n^i ^(^). (4.6) 

7 




Fig. 4.1. Typical trajectory 



The trajectory is sketched in Figure 14.11 

When the initial point x is in D2 or D3 we just need to consider -/V(l), such that (i?(l), Z{1)) £ 
D 1 and Z{2) = 7 2 z - r(l) - iV(l) e B!>(~f T ), for some T greater than a suitably chosen To, where 

( r + ( + \(l-j)z, for xeD 2 , 

r(l) = <^ r* + A(l - 7)2!, for x e D 3 , 

{ r - £ + \(l -y)z, for x e D 4 . 

The rest of the analysis follows. □ 

4.2. Proof of Lemma 13.31 Assume without loss of generality that C n D\ has positive 
Lebesgue measure. Let us show that the set C\ = C V\D\ satisfies the "smallness" property in 
the statement. 

In Lemma 13.21 we essentially proved that we can reach any bounded Lebesgue-measurablc 
set B of positive measure from any state 1 in a finite number of steps. However, in order to 
give an upper bound on the required number of steps, we defined the set E — I x [0, a] and we 
introduced the measure ifE, which is defined as the Lebesgue measure of the set obtained via 
intersection with E. We obtained that for any < e < 1, for any B such that <pe(B) > 0, and 
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for any x, there exists T$(B, e,x) > 0, such that for any T > To(_B,e,x), there exists a strictly 
positive constant K = K(T,B,x) > such that 



F T (x,B) >K{T,B,x)e<p E (B) 



(4.7) 



In order to prove smallness, we need to eliminate the dependence on the specific starting 
point x and on the destination B of K and To. 

The dependence on x can easily be dealt with. Fix some destination set B. Since the set C\ 
is bounded, we can prove that both 

r (B,e,Ci) := sup T {B,e,x) and K(T, B, d) := inf K(T,B,x). 



arc finite. Indeed, revisiting the proof of Lemma 13.21 we have 



T (B,e,d) = sup 

ieCi 



1 



lo 



log 7 °z(l) 



1 



lo 



log 7 6 + 7|inf J| 



< oo. 



(4.8) 



Furthermore, in 



we have that 



p 2 (T,B,e,C 1 ) := inf p 2 {T,B,e,x) > 
xeCi 



by the fact that C\ is bounded. For 2 < t < T + 1, the p t > do not depend on the starting 
point. Again, via the boundedness of Ci, we have that there exists a constant 77 > such that 
for all x S Ci, we have h x (v, w) >-q for all (v, w) £ B s (-/ T ). 

The dependence on the destination set B is more tricky. The highest inconvenience is the 
dependence on B of the proposed value for To, which, as seen in (|4.8[) . depends logarithmically 
on S. Hence, for any e, one can always find a set B which is such that S needs to be arbitrarily 
close to in order to satisfy (|4.2[) . This in turn leads to To being arbitrarily large (whereas we 
would like it to be constant). This is due to the deterministic geometric convergence of Z(t) 
towards in region D3. In other words, if starting from a strictly positive value, Z(t) cannot 
reach in a finite time following the trajectory we considered. 

Let us instead pick a finite closed interval / and constants A > <5 > and let us define the 
set A := I x [6, A] . Clearly, v{A & ) = v(A), so in this case we are allowed to pick e = 1. We also 
define r)(T, A, C\) := va-i x& c u {v,w)eA{i T ) h x {v,w) > and 



P2 (T) := inf 

re[7z(l)-A7 T . 7 2(l)-<57 T ] 



r*+<>-C-(l+A)r-A(l-7)z(l) 



N{dn 2 ). 



r*-2f-(l+A)r-A(l-7)z(l) 



With T = T (A, 1,Ci) given by (|4.8[) , but for the chosen 5 > (no dependence on any 
specific destination set), for all T > Tq we define 



nT+l 
t=2 Pt 



PA, 



where ipA is defined in (|4.1[) . By (|4.6[) we can conclude that C\ is i/^-small. 
In a similar fashion, we can prove that the entire set C is small. □ 
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4.3. Proof of Theorem 13.11 In the following we define a Lyapunov function V which is 
unbounded off petite sets, that is for any n < oo the sublevel set Cy(n) := {y : V(y) < n} is 
small. We show that there exist constants a, b, c > and a set C = [—a, a] x [0, b] such that the 
drift satisfies 

VV{x) := E X V(X(1)) - V(x) < -1 + ct c . 

A consequence of this result is that, by Theorem 9.1.8 of [§] and by Corollary 13.51 which shows 
that the set C is small, the chain is Harris recurent. Furthermore, by Theorem 10.0.1 of [5], it 
admits a unique invariant measure it. Finally, by Theorem 13.0.1 of [9] and by Corollary 13.41 we 
get the finiteness of it and we conclude. 

Let us now delve into the details: 

We prove that the function JJ:Ix H + — > R + , 



i ^ H(x) = (r + Xz) 2 + (r + (A + y)zf 



(4.9) 



is a Lyapunov function for the system (|2. 14112. 15|) . 

Since /j ^ 0, matrix A\ can be diagonalized as A\ = M±A\M^~ , where: 



Mi 



-A — /i 
1 



-A 



Mi" 



-1 -A 
1 A + /x 



, Ax 



1 

1 - ,u 



Then, for G £>i, we can write 

X(i + 1) = AI 1 A 1 M^ 1 X(t) + Co + JV (* + 1). 
Via the change of variable Y(t) = X{t), we have that 

Y(t + 1) = A 1 Y(t) + Mf 1 (Co + iV (t + 1)). 

Denote 



n*) = 

and Ai := {i/ 6 R 2 : Mxj/ £ Z?i} or, more explicitly. 



U(t) 
V(t) 



Ai = {(«, u) : m + v > 0, (A + + Xv > 0}. 
Then, for Y(t) G Ai, the system can be rewritten as: 

U(t + l) = U(t)-±(N(t + l) + 0, 
V(t + 1) = (1 - fx)V(t) + I(JV(t + 1) + C). 



(4.10) 



(4.11) 



Hence, for a point x = (r,z) G D\, the Lyapunov function can be rewritten in the new 
variable y = M7 x x = (u, v) G Ai, u = — — (r + Xz) and v = —(r + (A + fj,)z) as 

Wi : Ai R+, := H(M iy ) = ^{u 2 + v 2 ). (4.12) 

Let us compute the drift of H in x G Z?i, T>H(x) = T>Wi(y): 

VW x {y) = E y (Wi(F(l))) - Wi(y) = -M 3 (2 - m)« 2 + 2C/i(l - fi)v - 2(fm + 2(C 2 + a 2 ). 
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Since [i > 0, the set Si of y G Ai for which VWi(y) < —1 is given by the hypergraph of the 
second degree polynomial function 



9i(v) 



M 2 (2 - AO 



v 2 + (1 - n)v + 



2(C 2 + o 2 ) + 1 



2C v r/ 2Cm 

namely Si = hyp + (<7i) n Ai = {(u,v) G Ai : u > gi(u)}. The complementary set Si on which 



VW\{y) > —1 is Si = hyp (gi) n Ai. Wc represent both sets in Figure 4.2(a) In the original 
coordinates, we denote the corresponding set by C\ = M\S\ = {M\y : y 6 Si}. 



v = 
















V 










(a) 5i 






r'-C 



(b) S 2 




r" +{ 



(c) 54 (d) C1UC2UC3U C 4 

Fig. 4.2. Positive drift sets 

Next, consider the case X(i) s £*2- The eigenvalues of A2 are {1, 1 — A — fi}. Let us consider 
the diagonalization A2 = M2A2M2 , where 



M 2 = 



1 -A 
1 



M2 1 = 



1 A 
1 



A 2 = 



1 

1 - X-/JL 



Let us apply the coordinate change Y(t) = M 2 X(t). For Y(t) G A 2 = {y : M%y G D2} we get 
y(i + 1) = A 2 Y(t) + N Q (t + 1) + Co- Explicitly, the domain is 

A 2 = {(u, u) : < u - \v < r* - (, v>0}. 

Equivalently, denoting U(t) = R(t) + XZ(t) and V(t) = Z(t), for (U(t), V(t)) G A 2 , we have 



U(t + 1) = U(t) + N(t + 1) + C, 
7(* + l) = (l-A-|i)V(t). 



(4.13) 
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We write the Lyapunov function (|4.9I) in the new coordinates: 

W 2 : A 2 -> R + , W a (y) := if(M 2 y) = u 2 + (u + ^vf. 

Then, the drift of H in a point x € D 2 can be written as VH(x) = T>W 2 (y), where y = (u, v) = 
M~ x x. Thus, 

VW 2 {y) = E y (W 2 (Y(l))) - W 2 (y) = 2(u + C 2 + a 2 - ^{2 - (A + M ))(A + /i)i' 2 
-2/i(A + p)ui> + 2Cit + 2/1(1 - A - [i)(v + C 2 + a 2 . 

Since Xv < u < Xv + r* — £, we have that 

VW 2 {y) < 2a 2 - 2( 2 + 4Cr* + 2C[2A + //(l - A - p)]t> - + H?{2 - m)« 2 

Hence, there exists a w + such that T>W 2 (y) < —1 for u > v + . Denote by S 2 ~ {y ~ (u,v) S 
A2 : v < v + } (which we represent in Figure 4.2(b) I and the corresponding set in the original 
coordinates by C 2 = M 2 S 2 = {M 2 y ■ y £ 52}. 

Since A4 = A 2l let us consider the case X(t) £ D4 and the same coordinate change Y(t) = 
M^X{t). For Y(t) £ A 4 = {y : M 2 y £ D±} we get Y(t + 1) = A 2 Y(t) + N {t + 1) - 
Explicitly, the domain is 

A 4 = {(u, v) : r* + £ < u - Xv, v > 0}. 

Equivalents, denoting again U(t) = R(t) + XZ(t) and V(t) = Z(t), for (U(t),V(t)) £ A 4 , we 
have 

' U(t + l) = U(t) + N(t + l)-t, . . 

V(t + l) = (l-\-n)V(t). [ ' 

We write the Lyapunov function (|4.9[) in the new coordinates: 

Wi : A 4 -> R+, Wi(y) := H(M 2 y) = u 2 + (u + fiv) 2 . 

Then, the drift of if in a point x £ D4 can be written as VH(x) = VW4(y), where y = (u, v) = 
M 2 _1 x. Thus, 

VW A {y) = m y {W A {Y{l))) - W 4 (y) = -2£u + £ 2 + a 2 - ^{2 - (A + /x))(A + fi)v 2 
-2/i(A + fi)uv - 2£u - 2^(1 - A - + £ 2 + a 2 . 

Since u > Xv + r* + £, we have that 

2W 4 (y) < 2cr 2 + 2£ 2 - 4£u - (2fx(X + fx)r* + 2£(x)v - /Lt(A + ^i) 2 (2 - /x)u 2 . 

Then, we have that VWiijj) < — 1 for points ?/ = (u,v) in the hypergraph of the second degree 
concave polynomial function 

g 4 (v) = ^ [2a 2 + 2£ 2 + 1 - (2 M (A + M )r* + 2^)v - /i(A + /i) 2 (2 - n)v 2 } . 



In Figure 4.2(c) we represent the set S4 on which VWi(y) > —1, namely S4 = hyp - (34) n A4. 
In the original coordinates, we denote the corresponding set by C4 = M 2 Si = {M 2 y : y g 54}. 
Finally, write the drift of the Lyapunov function in a point x = (r, z) £ D3: 

VH(x) = E(Az + r* + N) 2 + E((l - fi)(X + fi)z + r* + N) 2 - (r + Xz) 2 - (r + (A + fi)z) 2 . 
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Wc have that 

VH(x) < 2(r* 2 +<t 2 - r 2 ) + 2r* (A + (A + - fi))z - (A + [i) 2 fi(2 - \x)z 2 . 

Thus, T>H(x) < — 1 holds for points x = (r, z) € D% found outside the ellipse defined by 

£ = {(r, z) : 2r 2 + a(z - /3a" 1 ) 2 = 1 + 2r* 2 + 2cr 2 + /3 2 a -1 }, 

where we denoted a = (A + n) 2 [i{2 — /i) and (3 = r*(A + (A + ^)(1 — /i)). We denote by 
C 3 =mt(£)nD 3 . 

We have found that T>H < — 1 outside the bounded compact set C\ U C2 U C3 U C4, which 
concludes the proof. □ 



4.4. Proof of Theorem 13.61 Notice that if [i < —A, then the Z coordinate of the Markov 
chain cannot decrease. Hence the chain is trivially unstable (it is not even ^-irreducible). 

In the case — A < fi < 0, we turn again to [5] for proving the claim. In this case the Markov 
chain (|2. 14112. 15]) is ^-irreducible. We need to find a function H which satisfies the hypothesis of 
Theorem 11.5.2 from [3] to show non-positivity. 

Like in the proof of Theorem l3.ll consider the change of coordinates Y(t) = M^ 1 X(t) : which 
yields evolution (|4~TT|) for Y(t) = (U(t),V(t)) £ A 1 defined in (|4~T0]) . The state space under 
these new coordinates becomes A = {(u, v) : u + v > 0}. 

Define H : A ->• R + , 
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v 



H(y) 



log(w) if v > 1, 
otherwise. 



(4.15) 



Let us show that H(Y) has finite increments in any point of the state space, namely that 

sup E y \H(Y(l)) - H(y)\ < +00. (4.16) 

Consider y = (u, v) 6 A such that v > 1. Then, since < A + /i < A, we have that necessarily 
y e Ai. We obtain: 



Ey\H(Y(l)) - H(y)\ 



i(v) 



log 1 - n 



C + n 

[IV 



N(dn) 



Af(dn) \ogv, 



where £i(v) = [i — ( — — fi)v. If A^(l) > £i(v), then by the definition of H we have that 
H(Y(1)) = 0, in which case we are left with the second term (since v > 1 implies logw > 0). 

For a given v, the argument of the logarithm in the integral of the first term is less than 1 
for n comprised between £o(v) := fi 2 v — ( and £\{v). Furthermore, it is lower bounded by -, and, 
since log a < a — 1 for a > 0, we can write 



E y \H(Y(l)) - H(y)\ < f° V log(l 

J —OO \ 



C + n 



Af(dn) 



< - 



Af(dn) log v + Af(dn) log v 
e (v) Ji^v) 

C r la ^ v ) n r°° 
— + / — JV(dn) + / Af(dn) log v 

M« J-OO Jtn(v) 



All the terms in the sum above are bounded for v > 1. The last term converges to as v goes 
to infinity, via two applications of the L'Hopital-Bcrnoulli rule. 
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Next consider y = (u,v) G A such that v < 1. Thus H(y) = 0. Wc distinguish 4 cases: 
Mi y G Di, M\y G D 2: M\y G D 3 , or M\y G ZA4. Equivalcntly, the 4 cases write as 



y G Ai 
2/ G A 2 
2/ G A 3 
2/ G A 4 



= {(u, u) : u + v > 0, (A + pi)u + \v > 0, v < 1}, 

= {(it, u) : it + v > 0, -r* + ( < (A + /Lt)u + Aw < 0}, 

= {(it, u) : u + u > 0, -r* - £ < (A + n)u + Xv < -r* + Q, 

= {(it, v) : u + v > 0, (A + At)u + Aw < -r* - £}. 



In the case y G Ai, we have 
E,|i/(F(1)) - #(y)| = f l[V) log[(l - + M _1 (C + n)]JV(dn) 

*/ — OO 




7V(dn)((l - + /i^C - 1) - a(2^)- 1 /V" 1 e" £ ^ 



Both terms are bounded and converge to as v goes to -co (the first term again by application 
of L'Hopital-Bernoulli, the second one trivially). 

In the case y G A2, it is straightforward to check that we necessarily have v < 0. Define 
£ 2 ( y ) : = n - C, + pt(A + n)(u + v) - fiv. We have that if JV(1) > i 2 {y), then H(Y(1)) = by 
definition of H(-). Now since it + v > and A + fi > 0, it must be that 

+ + < 0. 

Then, 

rMv) 

E V \H(Y{1)) - H(y)\ = / log[-(A + n)(u + v) +v + /T 1 ^ + n)]M{dn) 

J —OO 

rMv) 

*/ — OO 

= / AA(dn)( M - 1 C-l)-(T(2^)" 1 /V" 1 e"^, 

*/ — OO 

which is clearly bounded. 

The case y G A4 is closely related: Define ti(y) := n + £ + /i(A + /i)(it + u) — /iv. We have 
that if 2V(1) > £4(2/), then i?(F(l)) = by definition of #(•)• Then, just like above, 

E y \H(Y(l)) - H(y)\ < - / M(dn)(p-^ + 1) - a(2 7 r)- 1 / V'e"^" , 

«/ — 00 

which is again bounded. 

In the case y G A3, it again holds that v < 0. Define £3(2/) := /tt — £ — (A + /i)(l — + u). 
We have that if iV(l) > £ 3 (y), then if(F(l)) = by definition of H(-). Now since u + v > and 
A + (j, > 0, it must be that 

- /i)(A + n)(u + v) + v < 0. 
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Then, 



E y \H(Y(l))-H(y)\= [' log^- 1 [(l-ii)(X + fx)(u + v)+C + nW(dn) 
< / A^(dn)/i- 1 C-(r(2^)- 1 / 2 ^ 1 e-^, 



which is again bounded. 

Let us now examine the drift VH(y) = W, y (H(Y(l)) — H(y) in some point y 6 Ai n {v > 1}, 
which wc denote by T>H(v) by abuse of notation: 

Mv) / C + n\ f°° 

VH(v) = log 1 - u + ^— M(dn) - / Mdn) log v 

J-oo V M« / Ji^v) 



AA(dn) log(l - M) + / log 1+ ^ r ) N(dn) 



N(dn) \ogv. 

ti(y) 



Taking the limit as v goes to infinity, we find that lim^^oo T>H{v) = log(l — fi) > 0. Indeed, 
each of the last two terms converge to 0. In other words, for any e > 0, there exists a vo(e) > 1 
which is such that for all v > vo, \T>H(v) — log(l — fi)\ < e. Taking a sufficiently small e, say 
eo = | log(l — yu), we get that 

W > v (eo), VH{v) > 0. (4.17) 

Define the set C :— {y — (u,v) : v < «o(eo)} H A. Then any y in the complementary 
C := C c n A is such that > sup vieC H (j/i) and Dif(j/) > 0. Using (|4TT5|) and (|4TTT)) we 

can apply Theorem 11.5.2 from [3] and conclude. □ 
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Appendix A. Impact of Delayed Heating. 

A.l. With Constant Coefficient of Performance. Assume the appliance is a heating 
system, which is used to maintain the temperature of a building at temperature T(t). We use 
the model in [51 Chapter III-E]. Let e be the coefficient of performance of the appliance, T(t) the 
room temperature, and 9(t) the outside temperature. Assume in this section that the coefficient 
of performance is constant. Assume the heating appliance consumes an amount of energy energy 
equal to d(t) at time slot t; the obtained room temperature T(t) is defined by the following 
equation 

d{t)e = K(T(t) - 6{t)) + C(T(t) - T(t - 1)) (A.l) 

where K is the building's leakiness and C its thermal inertia. Note that this equation is valid if 
wc assume that the appliance is used for heating only and not for cooling. We make no particular 
assumption on d(t), which typically depends on thermostat regulation by the building occupants. 

Consider two scenarios. In the first scenario, there is always enough supply and the appliance 
does not need to adapt. The energy consumed by the appliance is D a (t), the naturally occurring 
demand. Let T*(t) be the temperature obtained. Thus 

D a (t)e = K(T*(t)-6(t)) + C(T*(t)-T*(t-l)) t = l...r (A.2) 

In the second scenario, energy supply is not sufficient at times 1, . . . , r — 1 and is sufficient at 
time r. The energy consumed at times 1, . . . ,r — 1 is D a (t) — F(t) where F(t) is the frustrated 
demand. The accumulated frustrated demand at time r — 1 is 

T-1 

Z{r-l) = Y J F(t) (A.3) 

i=l 

The obtained temperature is T(t) < T*(t) for t = 1, . . . , r - 1. We assume that T(0) = T*(0). 
Assume that at time r, we use the energy required to obtain the same temperature as in the 
first scenario; let D a (r) + Z(t) be this amount of energy, i.e. Z(t) is the amount of backlogged 
demand that remains at time r. We have thus 

(D a (t) - F(t))e = K(T(t) - 0(f)) + C(T(t) -T(t-1)) t = 1 . . .r - 1 (A.4) 

(D"(t) + Z{r))e = K(T*(t) 6(t)) + C(T*(r) - T(r 1)) (A.5) 

By summation of the last two equalities it comes: 



D a {t) + Z( T )-Z{T-1) 



-l 



= K y^{T{t) 9{t)) + (T*(t) - 9(r))j + C (T*(r) - T*(0)) (A.6) 
and by summation of Eq. (j A.2|) , 

($> a (t) ) e = K ( !>*(*) - <?(*))) + C (T*(r) - T*(0)) (A.7) 



\t=i 

thus, by subtraction: 



Z(r) Z(r - 1) = -f ( ]>>*(*) - T(t)) ) < (A.. 



Therefore, for this simple model, the backlogged demand diminishes, i.e. the total energy con- 
sumption is reduced if heating is delayed. In other words, there is positive evaporation. 
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A. 2. With Heat Pumps. The model above assumes that the coefficient of performance e is 
constant. For some systems, this does not hold as it might depend on the initial room temperature 
T(t), the target room temperature T*(t) and the outside temperature 0(t). In particular, for 
heat pumps, the coefficient of performance is reduced when the difFcrcnccs T*(t) — T(t) and 
T*(t) — 6(t) are large. To understand what happens with such systems, we revisit the second 
scenario above and, for simplicity, consider the simpler case where all demand is frustrated at 
times 1 to r — 1. Equations (|A.4[) and (|A.5[) are now replaced by 

= K(T(t) - 6(t)) + C(T(t) - T(t - 1)) < = l...r-l (A.9) 
(D a { T ) + Z{r))e'{T) = K(T*(t) - 6(t)) + C(T*{t) - T(t - 1)) (A.10) 

and D a (t) = F(t) for t = 1, . . . ,r— 1. Here we assume for simplicity that e is maintained constant 
in scenario 1 but the coefficient of performance is smaller for scenario 2, i.e. e'(r) < e. This is 
because the heat pump needs to produce warmer water at time r than in scenario 1 to be able to 
reach the target room temperature, and if the outside temperature is low, the efficiency is less. 
Summing these equations and combining with Eq. (|A.7[) gives 

Z(r) - Z(t " 1) = ~f (£V W - T(t))\ +(l- ^) 0D» + Z(r)) (A.ll) 



Compare with Eq. (|A.8[) : the last term is positive, and if large, it may well be that Z(j) > 
Z(t — 1). i.e. there may be some negative evaporation. 



